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ABSTRACT 

We examine the clustering properties of a population of quasars drawn from fully 
hydrodynamic cosmological simulations that directly follow black hole growth. We find 
that the black hole correlation function is best described by two distinct components: 
contributions from BH pairs occupying the same dark matter halo ('1-halo term', 
CBH.ih) which dominate at scales below ~ 300kpch~^, and contributions from BHs 
occupying separate halos ('2-halo term', ^BH,2h ) which dominate at larger scales. 
From the 2-halo BH term we find a typical host halo mass for faint-end quasars (those 
probed in our simulation volumes) ranging from M ~ 10^^ to a few IO^^Mq from z = 5 
to 2 = 1 respectively (consistent with the mean halo host mass). The BH correlation 
function shows a luminosity dependence as a function of redshift, though weak enough 
to be consistent with observational constraints. At small scales, the high resolution of 
our simulations allows us to probe the 1-halo clustering in detail, finding that ^sH.ih 
follows an approximate power law, lacking the characteristic decrease in slope at small 
scales found in 1-halo terms for galaxies and dark matter. We show that this difference 
is a direct result of a boost in the small-scale quasar bias caused by galaxies hosting 
multiple quasars (1-subhalo term) following a merger event, typically between a large 
central subgroup and a smaller, satellite subgroup hosting a relatively small black 
hole. We show that our predicted small-scale excess caused by such mergers is in 
good agreement with both the slope and amplitude indicated by recent small-scale 
measurements. Finally, we note the excess to be a strong function of halo mass, such 
that the observed excess is well matched by the multiple black holes of intermediate 
mass (10'^ - 10^ Mq) found in hosts of 7\f - 4 - 8 x 10" Af©, a range weU probed by 
our simulations. 



1 INTRODUCTION 



With supermassive black holes beinR found at t he cen- 
tre of most galaxies (|Korniendv fc Richstonelll995h . inter- 
est in quasars has increased significantly, with substantial 
investigation into fundamental relations b etween black hole 
mass es and their host galax i es' pro p erties (IMagorrian et al.l 



1998!; 'Ferrarese & Merritt' '2000'; 'Gebhardt et al.l l2000l : 
Trcmainc ct al. 2002; Graham & Driver 2007). In addition 
to these relations, statistical studies of the spatial cluster- 
ing of quasars provide the potential to better understand 
the relation between quasars, their hosts and the underly- 
ing dark matter distribution, as well as estimate quasar life- 
times (see, e.g., [Haiman & Hui 2001; Martini & Wein berg 
I2001I ) across a relatively large range of redshift. For exam- 
ple, strong clustering would suggest quasars should reside in 
massive groups. If so, they should be rare and in order to re- 
produce the quasar luminosity density, they must have long 
lifetimes. Conversely, low correlation would suggest more 
common quasars, and thus shorter quasar lifetimes. 

Early studies of quasar clustering produced varying re- 



sults for the clustering amplitude, with no clear agreement 
on overall evolution with redshift, som e suggesting mini- 
mal or decreasing clustering evolution (|Mo fc Fang Il993l : 
ICroom fc Shankd |l996|), w hile others found an increase 
in cl ustering with redshift (|Kundia 119971 : iLa Franca et al.l 
Il998l ). These findings were generally poorly constrained due 
to the small sizes of available quasar samples. With the 
emergen ce of large scale surveys such as Sloan Digital Sky 
Survey (lYork et all i2000) and the Two-degree Field QSO 
Redshift Survey ( Lewis et al.ll2002l ') , substantially larger cat- 
alogs have been compiled, permitting more detailed in- 
vestigation into the clustering properties of quasars, and 
many recent studies have been made into t his area (e.g. 



La Franca et al.l 



20051: Shen et all 



2008 



19981: iPorciani et all [20041; ICroom et all 



20071: iMyers et al.ll2q07al: Ida Angela et all 



Shen et al.l l2009bl : IRoss et al.l l200a ). These recent 



studies have found evi dence for an increase in clustering am- 
plitu de with redshift (jLa Franca et al.lll998l : IPorciani et al.l 
[2OOJ), primarily for z > 2, in agreement w i th prediction s 
from simulations (see, e.g. lBonoli et al.ll2009l : ICrotonll2009l '). 
In addition to overall evolution, the luminosity depen- 
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dence (if any) of large-scale clustering can provide significant 
insight into what quasar populations domina te different lu- 
mino sity ranges. For example, the model of iHopkins et alJ 
l|2005 a.b.c d. 2006 ) suggests that bright and faint quasars 
are similar objects which are observed at different phases 
of their lifetimes, rather than being fundamentally differ- 
ent populations of quasars (as simpler, 'on-off' models as- 
sume). This model would suggest that both bright and faint 
quasars should populate similar halos. Thus, while there 
may be some correlation between peak luminosity and host 
halo mass, clustering dependence on instantaneous lumi- 
nosity should be relatively weak, particularly when com- 
pared to more tra ditional 'on-off' models of quasar lumi- 
nosity (jLidz et al.ll2006. ). Recent observational studies have 
generally found a lack of lu minosity dependen c e in the cor- 
relation function (see, e.g., Croom et alj 20051: iMvers et al. 
l2007al : Ida Angela et al.l |2008| ). though IShen et al. (2009b) 
found evidence for some, though weak, luminosity depen- 
dence. Several semi-analytic models have also been used, 
finding differing luminosity dependences, such as a signif- 
icant dependence for sufficient luminosity ranges, but lim- 
ited when considering o nly luminosities probed by observa- 
tion IjBonoli et al.ll2009l) , or weak depend ence at low re dshift 
{z < 1), but stronger at higher redshift (|Crotonll2009r ). 

In addition to large scale behavior, the possibility of 
excess quasar clustering on very small scales has arisen in 
several recent studies. While some observed quasar pairs are 
believed to be the result of gravitationally lensed quasars, 
it has been proposed that others may be physically distinct 
quasar binaries, which would suggest quasars cluster much 
more strongly on small scal es than extrapolation of large 
scale clustering would imply (lDiorgovskilll99ll : Hewett et al.l 
1 19981 : iKochanek eraLlll999l : iMortlock et al.lll999l) . suggest- 
ing a connecti on between g alaxy mergers and quasar ac- 
tivity (see, e.g. iKochanek et al. 1999). However, investigat- 
ing the smallest scale clustering has typically been prob- 
lematic due to observational limitations (such as fiber col- 
lisions preventing small-separation pairs from being distin- 
guished as distinct objects) and sample sizes insufficient for 
probing the smallest scales, where quasar pairs are rare. 
There have been several studies probing clustering at sub- 
Mpc scales, generally finding no excess clustering relative to 
an extrapolation of the large-scale clustering beh avior (see, 
e.g. IShen et al. 20 09al : IPadmanabhan et al.ll2009l ). However, 
these studies have been limited to scales above 100 kpc h~^, 
while several recent studies have managed to probe even 
smal ler scales, wh ere they do indeed find a significant ex- 
cess l|Hennawi et a l. 2006; Myers ct al. 2007b, 2008). 

In particular, Henn awi et al.l (|2006l ) studied binary 
quasars from SDSS and 2dF Quasar Survey to compute 
quasar clustering for scales as small as 20 kpc h~^ (co- 
moving), and found significant excess clustering relative 
to the large scale extrapolation (by an order of magni- 
tude at comoving scales below 100 kpc h~^, and growing 
stronger with decreasing scale). This excess implies that the 
quasars are more strongly clustered than galaxies at these 
small scales, supporting the theory that quasar activity is 
trigg ered by galaxy inter a ctions. Using the qu asar sample 
from iMvers et al.l (|2007al ). iMvers et al.1 (|2007lj) found only 
a slight excess in small-scale clustering, and put an upper 
limit for the excess at a factor of 4.3+1.3 for physical scales 
of r^ 28 kpc h~^. They suggest that the significantly larger 



excess of iHennawi et al.l (|2006l ) is a result of a selection ef- 
fect, possibly due to studies tending to target tracers of the 
Lya forest, causi ng a bias tow a rd z > 2, which may be more 
highly clustered. iMvers et al.l (120081 ) used a complete spec- 
troscopic sample of quasars over physical scales of 23.7-29.9 
kpc h~^ from SDSS to find an exces s clustering f actor o f 
~4, consistent with the upper limit of lMvers et al.l (l2007a), 
which, while 2a below the excess found by iHennawi et al.l 
(2006 ). nonetheless supports the general finding of a clus- 
tering excess which may be a result of galaxy interactions. 

In this paper, we use cosmological hydrodynamic sim- 
ulations which directly model the growth, accretion, and 
feedback processes of black holes to investigate the proper- 
ties and underlying causes of black hole clustering. Although 
the simulation volume limits our analysis to black hole lumi- 
nosities and host group masses below those typically stud- 
ied, the self-consistent modeling of black holes allows us to 
study the clustering behavior without post-processing mod- 
els. Additionally, the high resolution allows us to investi- 
gate clustering behavior at extremely small scales, well be- 
low those studied with semi-analytic models, thereby provid- 
ing a means of using simulations to investigate the observed 
small-scale excess for the first time, and provide a physical 
explanation for the underlying cause. 

In Section 2 we describe the numerical modeling for 
the black holes formation and accretion (Section 2.1) the 
simulation parameters used (Section 2.2), the details of the 
subgroup finder (Section 2.3) and our method of calculating 
correlation functions (Section 2.4). In Section 3 we investi- 
gate the quasar clustering properties at both large and small 
scales, and we summarize our results in Section 4. 



2 METHOD 

2.1 Numerical simulation 

In this study, we analy se the set of simulations published in 
iDi Matteo et al.l (|2008l ). Here we present a brief summary 
of the si mulation code and the method used. We refer the 
reader to lDi Matteo et al.1 (|2008l ) for all details. 

The code we use is the massively parallel cosmological 
TreePM-SPH code Gadget2 (Springel 2005), with the ad- 
dition of a multi-phase modeling of the ISM, which allows 
treatment of star formation (Springel & Hernquist 2003), 
and black hole accretion and associated feedback processes 
(Springel et al. 2005, Di Matteo et al. 2005). 

Black holes are simulated with coUisionless parti- 
cles that are created in newly emerging and resolved 
groups/galajdes. To find these groups, a friends-of-friends 
group finder is called at regular intervals on the fly (the time 
intervals are equally spaced in log a, with A log a — log 1.25), 
finding groups based on particle separations below a speci- 
fied cutoff. Each of these groups that does not already con- 
tain a black hole is provided with one by turning its dens- 
est particle into a sink particle with a seed black hole of 
fixed mass, M = 5 x 10^ h~^ Mq. After insertion, the black 
hole particle grows in mass via accretion of surrounding gas 

according to Mbh = ,^2, 2^3'/ 2 (JHovle fc Lvttletonlll939l : 
iBondi fc Hwi3 Il944l : iBoridil Il953 'l. and by merging with 
other black holes. Note that within the simulations, it is 
assumed that accretion is limited to a majdmum of 3 times 



© 200? RAS, MNRAS 000, ??-?? 



Quasar Clustering in Cosmological Hydrodynamic Simulations 3 




Figure 1. An example of the distribution of black lioles in tlie simulations: The same slice (2 Mpc h~^ thick) through the D6 simulation 
at z=l,2,3,4. The positions of black holes in different luminosities bins (L < 10*1/0 - Orange; 1O*-L0 < L < IO^Lq - Pink; IO^Lq < 
L < 10^'' -Lq - Blue; L > 10^" Lq - Green.) arc plotted on top of the gas density distribution (shown in the the gray scale). 



the Eddington rate, although very few sources accrete above 

MEdd- 

The accretion rate of each black hole is used 
to compute the bolome tric luminosity, L = rjMBiic^ 
IShakura fc Sunvaevll973l '). Here rj is the radiative efficiency, 
and it is fixed at 0.1 throughout the simulation and this 
analysis. Some coupling between the hberated luminosity 
and the surrounding gas is expected: in the simulation 5 per 
cent of the luminosity is (isotropically) deposited as ther- 
mal energy in the l ocal black hole kernel, acting as a form 
of feedback energy (|Di Matteo et al.lbOOsI '). 



Table 1. Numerical Parameters 



Run 



Boxsizc 
h~^Mpc 



N^ 



m-DM 



/i-IMq h-^MQ h-ikpc 



D6 
E6 



33.75 
50 



2 X 486^ 
2 X 486^ 



2.75 X 10^ 
7.85 X 10'' 



4.24 X 10^ 
1.21 X lO'' 



Np : Total number of particles 

rrioM- Mass of dark matter particles 

mgas : Initial mass of gas particles 

e: Comoving gravitational softening length 



2.73 
4.12 



2.2 Simulation parameters 

Two simulation runs are analysed in this paper to allow for 
different volume size and resolution. The main parameters 
are listed in Table [T] Both runs were of moderate volume, 
with boxsizes of side length 33.75/i~^Mpc (D6 simulation), 
and 50/i"^Mpc (E6). For both simulations iVp = 2 x 486^ 
particles were used. The moderate boxsizes prevent the sim- 



ulations from being run below z ~ 1 to keep the fundamental 
mode linear, but provide a large enough scale to produce sta- 
tistically significant quasar populations. The limitation on 
the boxsizes is necessary to allow for appropriate resolution 
to carry out the subgrid physics in a converged regime (for 
further details on the s imulation methods, par ameters and 
convergence studies see lDi Matteo et al.l (|2008r )). 
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Figure 2. Relation between masses of dark matter halos and their most massive black holes. Color represents bolometric luminosity of 
the massive BH. 



2.3 Subgroup finder algorithm 

In addition to the on-the-fly friends-of-friends algorithm 
used to identify groups, a modified version of the SUB- 
FIND algorithm (|Springel et al-lbOOlT ) was run on the FoF- 
identified groups to determine the component subgroups 
(i.e. galaxies) within each group. These subgroups are de- 
fined as locally overdense, self-bound particle groups. To 
identify these regions, the algorithm sorts the particles 
within the parent group by density, and then analyzes each 
particle in order of decreasing density. For each particle i, 
the density of the 32 nearest neighbors are checked. If none 
are denser than particle i, it forms the basis for a new sub- 
group. If a single particle denser than i is found, or if the 
closest two denser particles belong to the same subgroup, 
particle i is assumed to be a member of that subgroup. If 
the two nearest particles denser than i are members of diflfer- 
ent subgroups, these two subgroups are stored as subgroup 
candidates, and are then joined into a new subgroup also 
containing i. After checking each particle in this manner, 
particles are checked for binding within their parent sub- 
group based on their position relative to the position of the 
most bound particle, and the velocity relative to the mean 
velocity of particles in the group. Any particle with positive 
total energy is considered unbound, and is removed from the 
subgroup, leaving the group divided up into its component 
subgroups (galaxies). 



2.4 Correlation Function 

To investigate the clustering properties of quasars, we use 
the two-point correlation function ^(r): 

dP ^ pl[l + i{r)]dV^dV2 (1) 

(|Peacocklll999l ). where dP is the probability of finding one 
object in each volume element dVi and dV2, separated by a 



distance r, with an average number density of po. We use the 
natural estimator ^(r) — -22 — 1 for computing ^, where DD 
and RR are the number of pairs of objects found with sepa- 
ration r in the simulation (DD) and in a random distribution 
of equal spatial density (RR) . For calculating RR, we used a 
random distribution of A'^^ = 6x 10^ objects to find the num- 
ber of pairs in a random sample, which is then normalized 
with a factor of (j^) (where No is the number of objects 
considered for the DD term) to correct for the increased 
spatial density of the random sources relative to the BHs in 
the DD term. Note th at the estimator ^(r) = dd-2dr+rr 
(jLandv fc Szalavlll993h has been shown to be more accurate 
(as it more effectively accounts for edge effects), but when 
considering small scales, both estimators provide equivalent 
results (iKer s cher e t al. 2000). Indeed, to confirm the valid- 
ity of the natural estimator, we compared results between 
the natural estimator and the Landy and Szalay estimator, 
and found that for the largest scales (> 5 h~^ Mpc) at low 
redshift, they differ by less than 5%, and the discrepancy is 
well below 1% everywhere else. 



3 RESULTS 

To illustrate the distribution of quasars (as a function of 
their luminosity) with respect to the underlying matter dis- 
tribution, in Figure [l] we plot a slice through the D6 sim- 
ulation at z = 1,2,3,4, with black hole positions indicated 
by colored dots for four luminosity range bins: L < 10* Lq - 
Orange; 10* Lq < L < lO'' Lq - Pink; IO^Lq <L< 1Q^° Lq 
- Blue; L > lO^^L© - Green. As expected, as supermassive 
black holes are hosted by galaxies, the quasars (particularly 
the most luminous sources) are located in some of the dens- 
est regions, with low redshift tending to exhibit more BHs, 
though at generally fainter luminosities. To characterize the 
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relation between black hole and host halo mass more pre- 
cisely, in Figure [2] we show the relation between the group 
halo mass and the mass of its most massive (central) black 
hole, with color representing the respective (instantaneous) 
quasar luminosity. There is a correlation between halo mass 
and BH mass, and to a lesser extent between halo mass and 
BH luminosity, with large halos tending to host more mas- 
sive, more luminous black holes than smaller halos, albeit 
with significant scatter. This is due to the lightcurve that a 
black hole has in our si mulations (regulated by the complex 
hydr odynamics, see e.g. lDi Matteo et al.ll2008l : iDegraf et al.l 
I2OI0I). We also find that as redshift decreases, the simulation 
is more densely populated with BHs, which tend to be more 
massive and less luminous than at earlier redshift. 

To study the relation between black holes and other 
structures, in Figure [3] we show the correlation functions 
of black holes found in the D6 (solid black) and the E6 
(solid pink) simulations for scales between lOkpch"^ and 
~ lOMpch"^ at z=l, 3, 5, with Poisson error bars. Note, 
the results from the two simulations are very similar, with 
the higher resolution D6 simulation showing a small boost at 
small scales (below ~ 200kpch~^). In general, we see £,bh 
typically takes the form of a power law (with some possible 
excess at small scales at 2: = 1). 

We also divide ^bh,d6 into two terms: a 1-halo term 
(dashed blue) produced by BH pairs occupying the same 
host group, and the 2-halo term (dashed green) produced 
by pairs occupying different groups. As expected, the 2-halo 
term dominates at large scales (above ~ 300 kpc h~^), while 
at smaller scales the 1-halo term dominates, indicating that 
our small scale clustering is really measuring BH properties 
within the scales of the host halos. A distinction between 
the 1-halo and 2-halo terms is expected (as BHs are hosted 
by galaxies within halos) a nd is consistent with t he theo- 
retical expectations (see, e.g. lCoorav fc Shethll2002l ). as well 
as what has been found in galaxy correlat ion functions (see, 
e.g. iMagliocchetti fc Porcianl.200a : .Zehavi et al.ll2004l 'l. 

3.1 Large Scale Clustering 

It may be expected that black holes will cluster similarly to 
their host galaxies (within their halos). To investigate the 
relation between BH clustering and that of their host ha- 
los, in the left column of Figure |4] we plot the 1-halo (blue) 
and 2-halo (red) contributions to the correlation function 
(at z = 1, 3, and 5) for BHs (solid lines) and galaxies (as 
identified by the subgroup finder described in Section 2.3) 
populating halos (i.e. groups) in the specified mass ranges 
(dashed lines). These mass ranges were chosen to reproduce 
the closest agreement between ^bh and ^subgroup in the 2- 
halo regime at each redshift, so as to be used as an indicator 
of the typical halo mass for BH hosts (at each redshift). The 
same is shown in the right column of Figure |4] where we 
only include some of the most luminous BHs in the simula- 
tions {10^ Lq < I/BH < 10^^1/0, a range which is probed by 
observations) . 

For the full BH population, the typical host mass in- 
creases slightly with decreasing redshift, from ~ W^^ Mq to 
somewhat below 1O^^M0 from 2 = 5 to 2: = 1 respectively. 
When limited to the luminosity range 10^ Lq < L < 10^" Lq, 
we again find increasing host mass with decreasing redshift, 
but with a sharper increase up to masses a few times 10^^ Mq 
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Figure 3. Two point correlation functions for the black holes in 
the D6 (solid black) and E6 (solid pink) simulations at z=l, 3, 5, 
with the 1-halo and 2-halo terms for the D6 simulation explicitly 
shown (dashed blue and green, respectively). 



at z = 1 (still in the faint end of the luminosity function, 
see iDegraf et all \2Qldi ) . 

We compare the typical host mass found in this way to 
the mean (median) mass of the host halos (see Table [2]) for 
several luminosity ranges, and find that the 2-halo clustering 
as described above does indeed provide an estimator for the 
mean host halo mass at the corresponding redshift in the 
simulations. In addition, the table shows that for a given 
halo mass the luminosity of its typical BH decreases with 
time, particularly at low redshift (below 2 ~ 2 — 3) , as seen 
more generally in Figure [2] This is shown explicitly in the 
bottom of Table [21 where we calculate the mean and median 
BH luminosities found within groups of specified halo mass 
ranges. Note that the mean quasar luminosity actually peaks 
at 2 = 3 for massive (M > 10^^ Mq) groups as a result of a 
few highly luminous sources. 

To better characterize the overall clustering strength, 
and in particular its luminosity dependence and evolution 
with redshift, we use the correlation length ro, defined as 
the scale at which ^(ro) — 1 [which we calculate using a lin- 
ear extrapolation of ^] . In Figure [5] we plot ro versus z for 
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Figure 4. Correlation functions for the D6 simulation BHs (solid) and subgroups within a specified mass range (dashed), at z=l, 3, 5, 
with 1-halo and 2-halo terms plotted separately (blue and red, respectively). The BH correlation function is plotted using all BHs (left) 
and using only those with 10® Lq < Lbh < 10^'' Lq (right). The mass range for ^group is chosen so as to find the closest agreement 
between ^bh 2h ^nd ^group. We also plot ^BH.ih using only BHs found in host groups in this fitted mass range (dot-dashed green line). 



BHs in several luminosity bins (solid colored lines) and, for 
comparison, groups in several mass bins (dashed grey lines). 
In general we find a weak evolution of the quasar clustering 
with redshift. This can be simply explained by the evolution 
of the bias of its underlying host halo masses. In particular, 
the correlation length for luminous (L > l(f Lq) BHs tends 
to decrease slightly as a function of decreasing redshift until 
2: ~ 3, following closely the bias of the 10^^ — 10^"^ Mq groups 
(consistent with the constraints on the host masses of these 
BHs). At fixed mass, these groups are less biased as a func- 
tion of decreasing redshift IJMo fc Whitell2002l : lBahcall et all 



[2OOJ). This is also in accord with our results from Figure |4l 
that the typical host halo mass remains roughly constant 
for z; > 3. For lower redshift (particularly z < 2), we instead 
see a significant upturn in ro versus z, corresponding to the 
increase in typical host halo mass, just as we found in Fig- 
ure |4] and Table m The lowest luminosity sources, however, 
show only minor change in ro , corresponding to a host mass 
which changes only slightly with redshift (consistent with 
the median host masses found in Table [2]). 

This luminosity dependence is sufficiently weak (less 
than a factor of 2 increase in ro across several orders of 
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Figure 5. Solid lines: Black hole correlation length as a function 
of redshift for several luminosity bins (colored lines) . Grey dashed 
lines : Group correlation length as a function of redshift for group 
mass ranges (from top to bottom) > W^^Mq, 10^^ — 10^'^Mq, 
5 X IQio - 5 X 10"Afo, IQio - 10" M^. 



magnitude in luminosity) to remain broadly consistent with 
the predictions from models that suRRe sting bright and 
faint quasars occu py similar halos (e.g. iLidz et al.l 120061 : 
iBonoli et al.l 12009 ). Indeed, our simulations produce com- 
plex lightcurves for our black holes, with luminosity vary- 
ing rapidly across several orders of mag nitude (see, e.g. 
iDi Matteo et al.l l2008l : lOeeraf et al.ll2010l ). This produces 
significant scatter in the relation between black hole lumi- 
nosity and host mass, so general agreement with lightcurve- 
based models is expected (which are indeed motivated by 
simulations similar to our own). Nonetheless, as seen in Fig- 
ure [2j there remains some correlation between BH instanta- 
neous luminosity and group mass, so a weak dependence on 
luminosity is expected even in this model. 



3.2 Small Scale Clustering 

Although we find that the 2-halo terms for BHs and sub- 
groups (galaxies) can be easily matched to provide a good 
estimator for typical host mass, there is significant discrep- 
ancy between their respective 1-halo terms (Figure [4] blue 
lines). The 1-halo BH correlation function is different both 
in shape and amplitude to the 1-halo term of galaxies, sug- 
gesting that, unlike at large scales, BHs do not cluster like 
their host galaxies on small scales. Or in other words, the 
distribution of BHs within halos does not follow closely that 
of their galaxies and hence does not trace the underlying 
matter distribution. 

In terms of amplitude, ^BH,ih can be adjusted by only 
considering the BHs in those groups that match the mass 
range constrained by the 2-halo term, thereby minimizing 
the suppression of ^BH.ih from the numerous BHs in groups 
too small to contribute to the 1-halo term (due to hosting 
only a single BH). As expected, in this case, (FigurelU green 
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Figure 6. Solid lines: The 1-halo BH correlation function at 
z=l,2,3 divided into components from BH pairs occupying sepa- 
rate subhalos (red) or co-habitating a single subhalo (blue), using 
the full population of BHs. Dotted lines: The 1-halo subgroup cor- 
relation function at z=l,2,3. 



line) the amplitude increases and is more in agreement with 
the 1-halo term of the subgroups (^subgroup, ih; at least at 
z=l-3 where the statistics are good enough). 

It is, however, hard to account for the substantial 
difference in shape: ^BH.ih follows an approximate power 
law, lacking the decrease in slope at small scales (below 
~ 200 — 300 kpc h~^) observed in ^subgroup, ih and expected 
fro m the 1-halo clustering produced by a general NFW pro- 
file llNavarro et al.lll996l : ICoorav fc Shethl2002l : IZehavi et all 
[2OOJ). Thus the BHs are distributed significantly differently 
than an NFW profile, showing a significant boost at small 
scales. 

We investigate the reason for this difference in the shape 
of the BH 1-halo term in terms of multiple BHs co-existing 
in a given subgroup. These BHs end up in a given subgroup 
as a result of mergers between their host galaxies, so that 
multiple BHs are expected to co-exist in a remnant (until dy- 
namical friction is able to bring them close enough together 
to eventually merge). 

To understand the effect this has on the small scale 
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for reference. 



clustering of BHs, we calculate the contributions to ^bh.iii 
from pairs of BHs occupying the same galaxy (we will call 
this the '1-subhalo' term) and from pairs of BHs occupying 
different galaxies within the same group ('2-subhalo' term), 
in analogy with dividing the overall correlation function into 
its 1-halo and 2-halo terms. We note that the existence of 
multiple BHs within a single subgroup necessarily indicates 
a previous merger event, since BH particles are not inserted 
into galaxies which already contain a BH particle, and thus 
any 1-subhalo contribution is inherently a result of previous 
galaxy mergers. 

In Figure [6] we plot the 1-subhalo (solid blue) and 2- 
subhalo (solid red) components of ^BH.ih, together with 
Csubgroup.ih for subgroups in groups within the mass ranges 
listed in Figure [4] (dotted line). The 1-subhalo term does 
indeed have a steeper slope than the 2-subhalo term, and is 
most significant at small scales. The 1-subhalo term is most 
dominant at low redshift, and by 2 = 1 it dominates the en- 
tire 1-halo term. This is a result of having increasingly large 



groups at low redshift which have also undergone a relatively 
large number of mergers. These indeed contain multiply- 
occupied subgroups (see Figure [7]). We also find that, if re- 
stricted to BHs within the same host mass range as the sub- 
groups, the 2-subhalo term of ^BH.ih matches ^subgroup, ih 
quite closely. Thus we find that, within sufficiently large 
groups (such that the simulation contains enough groups 
hosting multiple BHs to produce a well-defined 1-halo term), 
CsH.ih has two distinct components: one due to BH pairs 
which occupy separate galaxies, exhibiting good agreement 
with ^subgroup. ih; and a steeper one caused by BH pairs 
which co-occupy individual galaxies as a result of previous 
galaxy mergers, causing a boost in the small-scale ^BH,ih, 
particularly evident at low redshift, where typical groups 
are largest and have undergone significant merging. 

In Figure [7] we plot the relative mass of each multiply- 
occupied subgroup and its host group, with circles indicating 
central subgroups, and triangles showing satellite subgroups. 
We clearly see that these multiply-occupied subgroups tend 
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to be the primary (central) subgroup within a given group, 
typically containing ^ 65—70% of the total group's mass. We 
also color-code the datapoints to show the number of BHs 
within each subgroup, and see that the central subgroup of 
larger groups tends to contain more BHs. 

To investigate the masses of BHs which populate these 
multiply-occupied subgroups, in Figure[H]we plot the mass of 
the largest (primary) black hole within a given subgroup rel- 
ative to the masses of the other BHs in the same subgroup, 
color-coded to show the number of black holes within the 
subgroup. In only a few rare cases do we have more than 
one massive BH, while in the majority of cases we have, at 
most, a single massive BH with one or more smaller black 
holes, generally within an order of magnitude of the seed 
mass. This suggests that the majority of BHs in multiply- 
occupied subgroups come from relatively small satellite sub- 
groups (hosting correspondingly small black holes) which 
have fallen in and merged with the large, central subgroup, 
but do not grow substantially, instead remaining much less 
massive than the primary BH in the galaxy. Additionally, we 
observe that over time the fraction of BHs in the simulation 
located within these multiply-occupied subgroups increases 
from 2% at z = 5 to 15% ai z — 1, as typical groups get 
larger and have had more opportunity for satellite subgroups 
to merge with the central subgroup. This increase in typi- 
cal group mass causes an increase in both the number of 
multiply-occupied subgroups, as well as an increase in the 
typical number of black holes found within them (as seen in 
Figure [8|, which produces the increased importance of the 
1-subhalo term with decreasing redshift seen in Figure |S] 

We will compare the small scale clustering from the 
simulations to observations in Section 3.4. 



3.3 Quasar Bias 

To further characterize the clustering properties of BHs, we 
now consider the quasar bias as a function of scale and red- 
shift. The bias is obtained by taking the square root of the 
ratio between ^bh and the DM correlation function (shown 
as dotted lines in FigureUJ. Based on our results of the small 
scale clustering, we expect the quasars to be strongly biased 
with respect to the DM distribution at small scales, par- 
ticularly at high redshift. This general trend is clearly seen 
in Figure |4j where ^dm (dotted lines) increases with time 
(due to gravitational collapse), while ^bh tends to decrease 
slightly (seen more clearly in Figure [5} . More importantly, 
we see that the BH clustering bias relative to that of DM is 
strongly scale-dependent, with ^bh exhibiting a significant 
increase in clustering at small scales (below ~ 300 kpch~^) 
due to the strong 1-halo term, whereas ^dm shows only a 
slight increase at these small scales. 

In the top of Figure [§] we plot the scale- 
dependent BH bias and subgroup bias (defined as 

\/^bh/Cdm; \/Csubgroup/(^DM, respectively) found within the 
hosts of the best-fitting mass ranges found in Figure [l] Here 
we see that the subgroup bias levels off (as did ^subgroup in 
Figures |4] and [6]), but the 1-subhalo term causes the BH 
bias to continue increasing to the smallest scales probed 
in our simulation. To show this more clearly, the middle 
of Figure [9] shows the bias of BHs relative to the sub- 

for z=l-5. Within a given host 





^?BH/?subgroup 
2-6.10" M„ hosts 



groups ^^BH/Csubg 



100 1000 10000 
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Figure 9. Top: Solid lines: Black hole bias defined as 
V^bh/'^dm, using only BHs occupying halos in the best-fitting 
mass ranges specified in Figure |4] Dotted lines: Subgroup bias 
defined as -v/^subgroup/CDM, using only subgroups occupying ha- 
los in the best-fitting mass ranges specified in Figure |4] Middle: 
Bias of BHs relative to subgroups ( -v/feH/Csubgroup) occupying 
halos in the typical mass ranges found in Figure |4] Bottom: Bias 
of BHs relative to subgroups (-v/^bh /^subgroup) occupying halos 
of mass 2-6 X 10" Mq. 



mass range, the BHs cluster very similarly to the sub- 
groups (galaxies), except at the smallest scales (below ~ 
100 kpc h^^), where we again see the increased clustering 
caused by the multiply-occupied subgroups remaining from 
merger events, as discussed earlier. Although we note that 
this small-scale bias appears to be redshift dependent, it 
is actually a result of the evolution of the host mass be- 
ing considered. At higher redshifts, the typical host mass is 
smaller, and thus fewer will have undergone subgroup merg- 
ers producing multiply-occupied subgroups (as confirmed in 
Figure [T]!, thereby making the small-scale boost less appar- 
ent. When considering behavior for a fixed group mass (as 
shown in the bottom of Figure |9]), we see that the bias be- 
tween BH and subgroup clustering is redshift-independent, 
and consistently exhibits a strong small-scale boost from 
past subgroup mergers. 
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Table 2. Mean (median) halo mass of parent group and mean (median) luminosity of daughter BHs in D6 simulation. 



z=5 



BH Luminosity 

All 

IO^Lq < Lbh < W'->Lq 



W^Lq < Lbh < 10"'L 







Mean (Median) Group Mass [10^°Mq] 

39.9 (10.1) 24.2(9.07) 19.3 (9.51) 14.8 (8.34) 12.6 (8.34) 

61.5 (18.7) 27.2 (9.66) 18.9 (8.80) 13.1 (7.53) 11.5 (7.12) 

252 (94.1) 54.9 (20.4) 38.2 (21.1) 25.5 (14.9) 15.9 (10.3) 



Group Mass 

Mgroup < 10"Mo 

1O"M0 < Afgroup < 1O"-5M0 

1O"-5M0 < Mgroup < lOl^Mg 

IQI^Mq < Mgroup < 10^^-^Mq 

1012-5 M© < Mgroup 



7.88 (7.69) 
8.66 (7.86) 
9.09 (8.23) 
9.45 (8.85) 
10.49 (9.51) 



Log(Mean (Median) BH Luminosity) [log(L(. 



8.96 (8.55) 
9.09 (8.76) 
9.64 (9.19) 
10.19 (9.81) 
10.33 (10.20) 
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Figure 10. Left: The projected correlation function from the D6 simulation, averaged across redshifts 1-3, and across 3 projected 
directions for the full BH population (solid line), for BHs found within groups of the typical host mass shown in Figure 2] (dashed black 
line), for BHs found in groups of mass 4 — 8 X 10^^ Mq (dashed pink line) , and for subgr o ups fo und in groups of mass 4 — 8 X 10^^ Mq 
(dotted pink line). We a lso plot the extensio n of the power law found in iPorciani et al.l (|200J) (step function) , and the observational 
results of lHennawi et al.l IJ2006) (asterisks) and lMvers et al.l ||2008| ) (triangle). Right: Same as left plot, but with H^p(iJmin, ^max) plotted 
for several lower-limits on the host group mass. 



3.4 Comparison with observations: Projected 
Correlation Function 



In order to compare with the observational constraints 
on the small sca le clustering (see Hcnnawi ct al. 20061: 
iMvers et al.ll2008l ). we compute the volume- averaged pro- 
jected correlation function Wp{Rniin, i?max). This projected 
correlation function is computed using the same estimator 
described in Section 2.4, but the separation between points 
is the projected separation onto the xy, xz, or yz plane, 
rather than the separation in three-space. Although these 
three projections provide comparable results, we average 
across the three directions to avoid any potential directional 
bias. We then average across redshifts 1-3 (to match the 



observational data redshift range), weighted by the num- 
ber of BHs at each redshift, an d plot the resul t in F igure 
1101 toge ther with the data from iHennawi et al.l (|2006l ) (as- 
terisks), |Mxers_et_alJ (|2008l ) (triangle), and the extension of 
the best-fit power law for the large scale clustering found by 
IPorciani et al.l (|200J) (step-function). We have also plotted 
the projected correlation function for subhalos found in our 
simulations for several host mass ranges (dotted lines). 

Figure [10] shows a remarkable agreement between the 
small scale clustering of BHs from the simulations with the 
observations. In particular, when considering BHs within 
groups in the mass range of 4 — 8 x IO^^Mq the small scale 
boost (magenta line) matches the observed clustering very 
well. For completeness we also show (dashed black line) the 
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signal expected from BHs in the hosts of the typical mass 
ranges shown in Figure |4] which is also in good agreement, 
although slightly lower normalization. Indeed, the observed 
small-scale excess can be explained as resulting from the 
merger-based boost found in our simulations, further em- 
phasizing the importance of such mergers on quasar evolu- 
tion. 

We also note that if the full BH population from our 
simulation is used, rather than those in the restricted mass 
range, we lose the small scale excess (solid line), since the 
majority of our BHs are found in groups too small to exhibit 
significant effects of subgroup mergers. 

To investigate the dependence of the projected corre- 
lation function on the host mass in more detail, in the 
right of Figure [TOl we plot Wp (i?min, flmax) for BHs hosted 
by groups with several different lower-mass cutoffs (from 
1 to 16 X 10^ ^Mq; colored lines), together with the observa- 
tional data. Here we see that including less massive groups 
causes an overall decrease in amplitude (as expected), and 
also suppresses the small scale excess, as a result of smaller 
groups being less likely to host a multiply-occupied sub- 
group. This suggests that, given sufficient observational 
data, small-scale clustering may provide a sensitive means of 
probing the typical mass of merging pairs of galaxies hosting 
supermassive black holes. As shown, the curves with a lower 
mass cut of -^ 4 — 8 x IO^'^Mq produce the best agreement 
with observation, implying that observed quasar pairs are 
typically located within groups of moderate size, compara- 
ble to those found within our simulation (and thus below 
the larger host masses typically associated with observed 
large-scale clustering) . 



4 CONCLUSIONS 

In this paper we have investigated the clustering of black 
holes within hydrodynamic cosmological simulations, its 
redshift evolution, luminosity dependence, and particularly 
the small-scale behavior. 

We have shown that the large scale clustering of black 
holes traces that of the galaxies within their host groups, 
and provides a predictor of the typical host mass, which 
for our simulations is found to be on the order of a few 
10^^ Mq. Although well below the typically found masses 
of ~ 2 X 10^^ - 10^^ Mq (iLidz et aLlbood : IRoss et al.|[2009l : 
iBonoli et al.ll2009l : [Shen et al.ll2009^ . this is consistent with 
our limited simulation volumes which can only follow the 
growth of the faint-end of quasar population (DeGraf et al. 
2010), and cannot follow formation of such massive groups. 
The typical host group mass shows some evolution with red- 
shift, most significant below 2; ~ 3, where typical host masses 
increase by up to a factor 10 (at z = 1). This low-redshift 
increase is distinctly luminosity dependent, with the more 
luminous sources (Lbh > IO^Lq) undergoing the most sub- 
stantial increase in typical host mass. Overall the evolution 
of clustering with redshift and luminosity is minor and con- 
sistent with current observational constraints (albeit in low 
luminosity populations this is yet to be fully constrained). 
The relatively weak dependence found in our simulations is 
consistent with the complex lightcurves we derive from our 
direct modeling in which quasar luminosities vary over rel- 
atively short timescales for a given host (as regulated by 



hydrodyn amical processes) . This is also consistent with the 
models of iLidz et all l|200d ). 

In addition to the large-scale clustering (the 2-halo 
regime), our simulations allow us to study the small scale 
clustering (the 1-halo term) of ^bh. We found that ^BH.ih 
follows a power law behavior all the way to the smallest 
scales. The clustering of black holes at small scale is unlike 
that of galaxies (or dark matter) . We showed that the 1-halo 
BH term can be subdivided into two components: 1-subhalo 
and 2-subhalo. The 1-subhalo term, fsubgroup.ih, represents 
the clustering of BHs within a galaxy and 2-subhalo that of 
BHs occupying different galaxies. We have shown that the 
1-subhalo is the one that provides the power law behavior, 
indicating that galaxies do contain multiple black holes as 
a result of mergers. These galaxies tend to be the central 
galaxy within relatively large groups (for our simulation), 
generally hosting at most a single massive BH with one or 
more smaller BHs, likely as a result of smaller satellite galax- 
ies merging with the large, central galaxy within the group. 
In the absence of these multiply-occupied galaxies, ^BH,ih 
and ^subgroup, ih exhibit very close agreement, but the inclu- 
sion of these merger remnants causes a significant boost in 
the small-scale BH clustering. This merger-based boost is 
most significant at low redshift, where typical group size is 
largest, though we find it in sufficiently massive groups at 
all redshifts. 

Though observational limitations make observing these 
scales difficult, several recent studies have found a small 
scale excess at scales below -^ 100 kpc h~^ IjHennawi et al.l 
[2OO6; Myers et al. 2008). The observed excess is in remark- 
able agreement to the one predicted by our simulations com- 
ing from groups approaching lO^^M©, which host mostly 
intermediate size black holes. This suggests that multiple 
black holes co-occupying a subgroup at low redshifts are 
likely faint (ish) AGNs hosted in Milky Way size halos that 
have recently undergone merging. We also note that galax- 
ies hosting multiple AGN IIKomossa et al.ll2003l: Gerke et al.l 
I2OO7I : iBarth et al.l l2008l : IComerford et al.l l2009bl) or inspi - 
ralling supermassive black holes ( Comerford et al.l [2009a|) 
have been found in recent studies, further supporting our 
conclusion of multiply-occupied subgroups. Although we 
leave more detailed investigation of the small scale BH pairs 
in our simulations (particularly with regard to the luminosi- 
ties of inspiralling black holes) for a future work, we note 
that our finding that multiply-occupied galaxies tend to host 
a single massive BH with one or more small BHs appears to 
be in keeping with the observation that most of the inspi- 
ralling BH pairs power only a single AGN (|Comerford et al.l 
'2009a). Given that, our agreement in small-scale merger- 
induced boost certainly reinforces the importance of galaxy 
mergers on the evolution of supermassive black holes. We 
also note this small-scale excess' sensitivity to the host mass 
suggests that future small-scale studies may provide a means 
to constrain the typical mass of merger events between 
galaxies hosting black holes, with current observational data 
combined with our simulations suggesting groups with typi- 
cal masses comparable to those probed in our simulations 
(from a few lO^M© to lO^^M©) produce the muhiply- 
occupied galaxies underlying the observed small scale ex- 
cess. 

We would like to point out however that there are sev- 
eral aspect of our modeling approach, including numerical 
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issues, in the simulations tiiat potentially affect our results 
on the small-scale clustering. We have a very simplistic pre- 
scription to determine how BHs merge with one another (im- 
posed by the limits on the resolution that can be achieved 
in these cosmological boxes) . The current prescription has a 
BH pair merge when BHs are separated by less than their 
smoothing length and if the BHs relative velocity is small 
compared to the local sound speed. Changes to this prescrip- 
tion could accelerate (postpone) BH mergers, which would 
result in a suppression (increase) of our small scale cluster- 
ing signals. It would be desirable to compare our results with 
other simulations which implement different prescriptions, 
or in the future to include more direct physical modeling of 
this region in higher resolution simulations. However, nei- 
ther of these are currently possible. A numerical issue that 
may affect the results of our one- halo term is that black holes 
need to be fixed to potential minima (calculated among the 
neighboring particles within the smoothing length used for 
the accretion model) in order to avoid them leaving their 
subhalo due to numerical N-body noise (and the fact that 
dynamical friction is hard to calculate for sink particles). 
However, in some instances this may cause a BH particle in 
a small subhalo in orbit in a bigger group to 'hop' to the 
potential minimum of the larger group. This effect may be 
exacerbated in situations where the small subhalo may be 
stripped of gas by infalling into a larger one. These effects 
could artificially increase the number of BHs within large, 
central halos, thereby boosting small scale clustering. How- 
ever, when we measure what fraction of BHs appear to 'hop' 
into the center of groups experiencing an unexpected jump 
in their position, we find that it is only ~ 1 — 2%. Future 
simulations and comparison amongst different approaches 
(once they become available) should of course attempt to 
characterize these effects more specifically. We further note 
however, that observational studie s have indeed found case s 
of galaxies hosting multiple BHs IjComerford et al.ll2009al ). 
so the existence of a one-subhalo term is expected. Addi- 
tionally, as seen in Figure 1101 the projected clustering of 
subgroups has a fundamentally different form than the ob- 
served quasar clustering. Thus the BHs cannot simply trace 
their host subgroups/galaxies and still produce the observed 
small scale excess, but rather a significant one-subhalo term 
is required to produce the small scale power law behavior. 

In future work we also plan to simulate larger volumes 
(which we are starting to be feasible with the most advanced 
technology) to allow us to study clustering of AGN at larger 
(mass and length) scales while simultaneously investigat- 
ing luminosity dependence for brighter sources more directly 
comparable to current and upcoming observational data, as 
well as providing increased statistics for the small scale clus- 
tering. 
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